%% it produces Figures 9, 10, 11, Table 4,

W_mat=welf_disc_all;

k_vec_base=k_vec;
W_mat_base=W_mat;
inde_mix = find(tau_vec==0);
[mm,inde_size_benefit]=max(W_mat(inde_mix,:));

W_opt = 1/cons_ss*W_mat(inde_mix,1:length(k_vec));
k_vec_opt=k_vec;

subsidy_opt.k_vec_opt = k_vec_opt;
subsidy_opt.W_opt = W_opt;
subsidy_opt.cons_ss=cons_ss;
str="/subsidy_opt.mat";
str_save=append(path_save,str);
save(str_save, 'subsidy_opt');

%% figures welfare
figure(100)
plot(100*(1-k_vec),100*1/cons_ss*(W_mat(inde_mix,1:length(k_vec)))-100*1/cons_ss*(W_mat(inde_mix,1)),'-b','Linewidth',6),hold on
%axis([0 50 -3 0.5])
xline(100*(1-k_vec(inde_size_benefit)),':k','Optimal size','Linewidth',2,'Interpreter','Latex','Fontsize',18,'LabelVerticalAlignment','bottom')
xlabel('$\vartheta+\tau$, \%','Fontsize',18,'Interpreter','Latex')
ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
grid on
le=legend('$\tilde\tau=0$');
set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthEast');
str1="/figures/figure_wefare_1.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)


figure(101)
plot(100*tau_vec(inde_solved(:,inde_size_benefit)==1),100*1/cons_ss*(W_mat(inde_solved(:,inde_size_benefit)==1,inde_size_benefit))-100*1/cons_ss*(W_mat(inde_mix,1)),'-b','Linewidth',6)
%plot(100*tau_vec,100*1/cons_ss*(W_mat(:,inde_size_benefit))-100*1/cons_ss*(W_mat(inde_mix,1)),'-b','Linewidth',6)
%axis([0 100 -1.5 0.5])
%vline(0,':b','Optimal Policy')
%xline(0,':b','Optimal','LineWidth',2)
yline(0,'--k','Status quo','Linewidth',2,'Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','left')
%xline(25,':r','ISS','LineWidth',2)
%yline(100*1/cons_ss*(W_mat(1,1)),'k--','Status Quo Welfare','LineWidth',2);
xlabel('$\tilde\tau=\tau/(\vartheta+\tau)$, \%','Fontsize',18,'Interpreter','Latex')
ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
grid on
le=legend('$\vartheta+\tau=0.35$');
set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthWest');
str1="/figures/figure_wefare_2.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)


figure(40)
plot(100*(tau_vec(inde_solved(:,inde_size_benefit)==1)),debt0_va_s_vec(inde_solved(:,inde_size_benefit)==1,inde_size_benefit),'-b','Linewidth',6),
%plot(100*tau_vec,debt0_va_s_vec(:,inde_size_benefit),'-b','Linewidth',6),
xlabel('$\tilde\tau$, \%','Fontsize',18,'Interpreter','Latex')
ylabel('Debt at entry: $b_0$','Fontsize',18,'Interpreter','Latex')
yline(debt0_va_s_vec(1,1),'--k','Status quo','Linewidth',2,'Interpreter','Latex','Fontsize',18,'LabelHorizontalAlignment','left')
le=legend('$\vartheta+\tau=0.35$');
set(le,'Interpreter','Latex','Fontsize',16,'Location','SouthEast');
settt
grid on
str1="/figures/figure_b0.pdf";
str_save1=append(root,str1);
print(gcf,'-dpdf',str_save1)


%% Tables at optimal policy
% benefits
Table1m=zeros(7,7);
Table1m(1,:) =  100/cons_ss*[welf_disc(inde_mix,inde_size_benefit)- welf_disc(inde_mix,1),NaN,welf_disc_s(inde_mix,inde_size_benefit)- welf_disc_s(inde_mix,1),NaN, welf_disc_n(inde_mix,inde_size_benefit)- welf_disc_n(inde_mix,1),NaN, welf_disc(inde_mix,inde_size_benefit)- welf_disc(inde_mix,1)-(welf_disc_s(inde_mix,inde_size_benefit)- welf_disc_s(inde_mix,1)+welf_disc_n(inde_mix,inde_size_benefit)- welf_disc_n(inde_mix,1))];
Table1m(2,:) =  100/cons_ss*[welf_disc_ss(inde_mix,inde_size_benefit)- welf_disc_ss(inde_mix,1),NaN,welf_disc_s_ss(inde_mix,inde_size_benefit)- welf_disc_s_ss(inde_mix,1),NaN,welf_disc_n_ss(inde_mix,inde_size_benefit)- welf_disc_n_ss(inde_mix,1),NaN,welf_disc_ss(inde_mix,inde_size_benefit)- welf_disc_ss(inde_mix,1)-(welf_disc_s_ss(inde_mix,inde_size_benefit)- welf_disc_s_ss(inde_mix,1)+welf_disc_n_ss(inde_mix,inde_size_benefit)- welf_disc_n_ss(inde_mix,1))];
Table1m(3,:) =  100/cons_ss*[y_disc(inde_mix,inde_size_benefit)- y_disc(inde_mix,1),NaN,y_disc_s(inde_mix,inde_size_benefit)- y_disc_s(inde_mix,1),NaN,y_disc_n(inde_mix,inde_size_benefit)- y_disc_n(inde_mix,1),NaN,  y_disc(inde_mix,inde_size_benefit)- y_disc(inde_mix,1)-(y_disc_s(inde_mix,inde_size_benefit)- y_disc_s(inde_mix,1)+y_disc_n(inde_mix,inde_size_benefit)- y_disc_n(inde_mix,1))];
Table1m(4,:) = -100/cons_ss*[i_disc(inde_mix,inde_size_benefit)- i_disc(inde_mix,1),NaN,i_disc_s(inde_mix,inde_size_benefit)- i_disc_s(inde_mix,1),NaN,i_disc_n(inde_mix,inde_size_benefit)- i_disc_n(inde_mix,1),NaN,   i_disc(inde_mix,inde_size_benefit)- i_disc(inde_mix,1)-(i_disc_s(inde_mix,inde_size_benefit)- i_disc_s(inde_mix,1)+i_disc_n(inde_mix,inde_size_benefit)- i_disc_n(inde_mix,1))];
Table1m(5,:) = -100/cons_ss*[conj_disc(inde_mix,inde_size_benefit)- conj_disc(inde_mix,1),NaN,conj_disc_s(inde_mix,inde_size_benefit)- conj_disc_s(inde_mix,1),NaN,conj_disc_n(inde_mix,inde_size_benefit)- conj_disc_n(inde_mix,1),NaN,  conj_disc(inde_mix,inde_size_benefit)- conj_disc(inde_mix,1)-(conj_disc_s(inde_mix,inde_size_benefit)- conj_disc_s(inde_mix,1)+conj_disc_n(inde_mix,inde_size_benefit)- conj_disc_n(inde_mix,1))];
Table1m(6,:) = -100/cons_ss*[chi_disc(inde_mix,inde_size_benefit)- chi_disc(inde_mix,1),NaN,chi_disc_s(inde_mix,inde_size_benefit)- chi_disc_s(inde_mix,1),NaN,chi_disc_n(inde_mix,inde_size_benefit)- chi_disc_n(inde_mix,1),NaN,  chi_disc(inde_mix,inde_size_benefit)- chi_disc(inde_mix,1)-(chi_disc_s(inde_mix,inde_size_benefit)- chi_disc_s(inde_mix,1)+chi_disc_n(inde_mix,inde_size_benefit)- chi_disc_n(inde_mix,1))];
Table1m(7,:) =  100/cons_ss*[sub_disc(inde_mix,inde_size_benefit)- sub_disc(inde_mix,1),NaN,sub_disc_s(inde_mix,inde_size_benefit)- sub_disc_s(inde_mix,1),NaN,-sub_disc_n(inde_mix,inde_size_benefit)+ sub_disc_n(inde_mix,1),NaN,  sub_disc(inde_mix,inde_size_benefit)- sub_disc(inde_mix,1)-(sub_disc_s(inde_mix,inde_size_benefit)- sub_disc_s(inde_mix,1)-sub_disc_n(inde_mix,inde_size_benefit)+ sub_disc_n(inde_mix,1))];


input.data =Table1m;
input.tableColLabels = {'Aggregate','','South','','North','','Reallocation'};
input.tableRowLabels = {'Total welfare','SS welfare only','Output','Investment','Congestions','Entrepreneur leisure','Net tax transfer'};
input.tableColumnAlignment = 'c';
input.dataFormat = {'%.2f'}; % three digits precision for first two columns, one digit for the last
input.tableBorders = 0;
input.tableCaption = 'Consumption equivalent gains at optimal subsidy: $\vartheta+\tau=0.35$, $\tilde{\tau}=0$';
input.tableLabel = 'tab:WelfareGains';
input.dataNanString = '';
input.makeCompleteLatexDocument = 0;
Table1 = latexTable(input);

str3="table1.tex";
str_save3=append(root,str3);
Table1_string=char(Table1);

writematrix(Table1_string,convertStringsToChars(str_save3),'FileType','text')

%% Tables at  other policy
inde_mix2 = find(tau_vec>=0.75,1);

% benefits
Table2m=zeros(7,7);
Table2m(1,:) =  100/cons_ss*[welf_disc(inde_mix2,inde_size_benefit)- welf_disc(inde_mix2,1),NaN,welf_disc_s(inde_mix2,inde_size_benefit)- welf_disc_s(inde_mix2,1),NaN, welf_disc_n(inde_mix2,inde_size_benefit)- welf_disc_n(inde_mix2,1),NaN, welf_disc(inde_mix2,inde_size_benefit)- welf_disc(inde_mix2,1)-(welf_disc_s(inde_mix2,inde_size_benefit)- welf_disc_s(inde_mix2,1)+welf_disc_n(inde_mix2,inde_size_benefit)- welf_disc_n(inde_mix2,1))];
Table2m(2,:) =  100/cons_ss*[welf_disc_ss(inde_mix2,inde_size_benefit)- welf_disc_ss(inde_mix2,1),NaN,welf_disc_s_ss(inde_mix2,inde_size_benefit)- welf_disc_s_ss(inde_mix2,1),NaN,welf_disc_n_ss(inde_mix2,inde_size_benefit)- welf_disc_n_ss(inde_mix2,1),NaN,welf_disc_ss(inde_mix2,inde_size_benefit)- welf_disc_ss(inde_mix2,1)-(welf_disc_s_ss(inde_mix2,inde_size_benefit)- welf_disc_s_ss(inde_mix2,1)+welf_disc_n_ss(inde_mix2,inde_size_benefit)- welf_disc_n_ss(inde_mix2,1))];
Table2m(3,:) =  100/cons_ss*[y_disc(inde_mix2,inde_size_benefit)- y_disc(inde_mix2,1),NaN,y_disc_s(inde_mix2,inde_size_benefit)- y_disc_s(inde_mix2,1),NaN,y_disc_n(inde_mix2,inde_size_benefit)- y_disc_n(inde_mix2,1),NaN,  y_disc(inde_mix2,inde_size_benefit)- y_disc(inde_mix2,1)-(y_disc_s(inde_mix2,inde_size_benefit)- y_disc_s(inde_mix2,1)+y_disc_n(inde_mix2,inde_size_benefit)- y_disc_n(inde_mix2,1))];
Table2m(4,:) = -100/cons_ss*[i_disc(inde_mix2,inde_size_benefit)- i_disc(inde_mix2,1),NaN,i_disc_s(inde_mix2,inde_size_benefit)- i_disc_s(inde_mix2,1),NaN,i_disc_n(inde_mix2,inde_size_benefit)- i_disc_n(inde_mix2,1),NaN,   i_disc(inde_mix2,inde_size_benefit)- i_disc(inde_mix2,1)-(i_disc_s(inde_mix2,inde_size_benefit)- i_disc_s(inde_mix2,1)+i_disc_n(inde_mix2,inde_size_benefit)- i_disc_n(inde_mix2,1))];
Table2m(5,:) = -100/cons_ss*[conj_disc(inde_mix2,inde_size_benefit)- conj_disc(inde_mix2,1),NaN,conj_disc_s(inde_mix2,inde_size_benefit)- conj_disc_s(inde_mix2,1),NaN,conj_disc_n(inde_mix2,inde_size_benefit)- conj_disc_n(inde_mix2,1),NaN,  conj_disc(inde_mix2,inde_size_benefit)- conj_disc(inde_mix2,1)-(conj_disc_s(inde_mix2,inde_size_benefit)- conj_disc_s(inde_mix2,1)+conj_disc_n(inde_mix2,inde_size_benefit)- conj_disc_n(inde_mix2,1))];
Table2m(6,:) = -100/cons_ss*[chi_disc(inde_mix2,inde_size_benefit)- chi_disc(inde_mix2,1),NaN,chi_disc_s(inde_mix2,inde_size_benefit)- chi_disc_s(inde_mix2,1),NaN,chi_disc_n(inde_mix2,inde_size_benefit)- chi_disc_n(inde_mix2,1),NaN,  chi_disc(inde_mix2,inde_size_benefit)- chi_disc(inde_mix2,1)-(chi_disc_s(inde_mix2,inde_size_benefit)- chi_disc_s(inde_mix2,1)+chi_disc_n(inde_mix2,inde_size_benefit)- chi_disc_n(inde_mix2,1))];
Table2m(7,:) =  100/cons_ss*[sub_disc(inde_mix2,inde_size_benefit)- sub_disc(inde_mix2,1),NaN,sub_disc_s(inde_mix2,inde_size_benefit)- sub_disc_s(inde_mix2,1),NaN,-sub_disc_n(inde_mix2,inde_size_benefit)+ sub_disc_n(inde_mix2,1),NaN,  sub_disc(inde_mix2,inde_size_benefit)- sub_disc(inde_mix2,1)-(sub_disc_s(inde_mix2,inde_size_benefit)- sub_disc_s(inde_mix2,1)-sub_disc_n(inde_mix2,inde_size_benefit)+ sub_disc_n(inde_mix2,1))];

input.data =Table2m;
input.tableColLabels = {'Aggregate','','South','','North','','Reallocation'};
input.tableRowLabels = {'Total welfare','SS welfare only','Output','Investment','Congestions','Entrepreneur leisure','Net tax transfer'};
input.tableColumnAlignment = 'c';
input.dataFormat = {'%.2f'}; % three digits precision for first two columns, one digit for the last
input.tableBorders = 0;
input.tableCaption = 'Consumption equivalent gains at optimal subsidy size: $\vartheta+\tau=0.35$, $\tilde{\tau}=0.75$';
input.tableLabel = 'tab:WelfareGainsISS';
input.dataNanString = '';
input.makeCompleteLatexDocument = 0;
Table2 = latexTable(input);
str3="table2.tex";
str_save3=append(root,str3);

Table2_string=char(Table2);
writematrix(Table2_string,convertStringsToChars(str_save3),'FileType','text')

%% IRFS at optimal policy
T_irf = 10/dt
Tm1=3;

figure(1)
plot(-Tm1+dt:dt:0,[100/cons_ss*W_time_s_vec(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
plot(-Tm1+dt:dt:0,[100/cons_ss*W_time_n_vec(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100/cons_ss*(W_time_s_vec(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100/cons_ss*(W_time_n_vec(1:T_irf,inde_mix,inde_size_benefit)),'-b','LineWidth',6),hold off
le=legend('South','North');
set(le,'Interpreter','Latex','Fontsize',16,'Location','southeast');
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
grid on
str1="/figures/figure_IRF_wefare.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)

figure(2)
% plot(-Tm1+dt:dt:0,[100/cons_ss*y_time_s_mat(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,[100/cons_ss*y_time_n_mat(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(y_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(y_time_n_mat(1:T_irf,inde_mix,inde_size_benefit)),'-b','LineWidth',6),hold off
plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'--r','LineWidth',10),hold on
plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100*log(y_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)/y_time_s_mat(end,inde_mix,1)),'--r','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100*log(y_time_n_mat(1:T_irf,inde_mix,inde_size_benefit)/y_time_n_mat(end,inde_mix,1)),'-b','LineWidth',6),hold off
%le=legend('South','Aggregate','North');
%le=legend('South','North');
%set(le,'Interpreter','Latex','Fontsize',16,'Location','southeast');
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Output, \% deviation from status quo','Fontsize',18,'Interpreter','Latex')
grid on
str1="/figures/figure_IRF_output.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)


figure(3)
% plot(-Tm1+dt:dt:0,[100/cons_ss*i_time_s_mat(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,[100/cons_ss*i_time_n_mat(end,inde_mix,1)*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(i_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(i_time_n_mat(1:T_irf,inde_mix,inde_size_benefit)),'-b','LineWidth',6),hold off
plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'--r','LineWidth',10),hold on
plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100*log(i_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)/i_time_s_mat(end,inde_mix,1)),'--r','LineWidth',6),hold on
plot(0:dt:T_irf-dt,100*log(i_time_n_mat(1:T_irf,inde_mix,inde_size_benefit)/i_time_n_mat(end,inde_mix,1)),'-b','LineWidth',6),hold off
%le=legend('South','North');
%set(le,'Interpreter','Latex','Fontsize',16,'Location','northeast');
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Investment, \% deviation from status quo','Fontsize',18,'Interpreter','Latex')
grid on
str1="/figures/figure_IRF_investment.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)

figure(4)
plot(-Tm1+dt:dt:0,[ones(1,Tm1/dt)]/dt,'--r','LineWidth',10),hold on
plot(-Tm1+dt:dt:0,[ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf-dt,(l_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
plot(0:dt:T_irf-dt,(l_time_n_mat(1:T_irf,inde_mix,inde_size_benefit)),'-b','LineWidth',6),hold off
%le=legend('South','North');
%set(le,'Interpreter','Latex','Fontsize',16,'Location','southwest');
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Mass of workers','Fontsize',18,'Interpreter','Latex')
grid on
str1="/figures/figure_IRF_labor.pdf";
str_save1=append(root,str1);
settt
print(gcf,'-dpdf',str_save1)

shift_exit_s=1;
shift_exit_n=1;
exit_ss1_s=exit_t_s(1,1,1);
exit_ss1_n=exit_t_n(1,1,1);

Tm1=3;
figure(5)
plot(-Tm1+dt:dt:0,100*shift_exit_s*[exit_ss1_s*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
plot(-Tm1+dt:dt:0,100*shift_exit_n*[exit_ss1_n*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_exit_s*([exit_ss1_s+exit_0_s(inde_mix,inde_size_benefit),exit_t_s(1:T_irf/dt,inde_mix,inde_size_benefit)']/dt),'--r','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_exit_n*[exit_ss1_n+exit_0_n(inde_mix,inde_size_benefit),exit_t_n(1:T_irf/dt,inde_mix,inde_size_benefit)']/dt,'-b','LineWidth',6),hold off
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Exit rate, in \%','Fontsize',18,'Interpreter','Latex')
%axis([-1 10 8.5 8.7])
%le=legend('South','North');
%set(le,'Interpreter','Latex','Fontsize',16,'Location','NorthEast');
str1="/figures/figure_IRF_exit.pdf";
str_save1=append(root,str1);
grid on
settt
print(gcf,'-dpdf',str_save1)

shift_entry_s=1;
shift_entry_n=1;
entry_ss1_s=entry_t_s(1,1,1);
entry_ss1_n=entry_t_n(1,1,1);

figure(6)
plot(-Tm1+dt:dt:0,100*shift_entry_s*[entry_ss1_s*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
plot(-Tm1+dt:dt:0,100*shift_entry_n*[entry_ss1_n*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_entry_s*([entry_t_s(1:(T_irf+dt)/dt,inde_mix,inde_size_benefit)']/dt),'--r','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_entry_n*[entry_t_n(1:(T_irf+dt)/dt,inde_mix,inde_size_benefit)']/dt,'-b','LineWidth',6),hold off
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Entry rate, in \%','Fontsize',18,'Interpreter','Latex')
%axis([-1 10 8.5 8.7])
%le=legend('South','North');
%set(le,'Interpreter','Latex','Fontsize',16,'Location','NorthEast');
str1="/figures/figure_IRF_entry.pdf";
str_save1=append(root,str1);
grid on
settt
print(gcf,'-dpdf',str_save1)

%% COMPARE SOUTH AT  at alternative policy with wrong timimg 
inde_mix2=find(tau_vec>=0.75,1);
T_irf = 10/dt
Tm1=3;
% 
% figure(7)
% plot(-Tm1+dt:dt:0,[100/cons_ss*W_time_s_vec(end,1,1)*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,[100/cons_ss*W_time_s_vec(end,1,1)*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(W_time_s_vec(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100/cons_ss*(W_time_s_vec(1:T_irf,inde_mix2,inde_size_benefit)),'-b','LineWidth',6),hold on
% le=legend('Optimal timing $\tilde\tau=0.00$','Suboptimal timing $\tilde\tau=0.75$');
% set(le,'Interpreter','Latex','Fontsize',20,'Location','southeast');
% xlabel('Year','Fontsize',18,'Interpreter','Latex')
% ylabel('Consumption gains, in \% dev. from  status quo','Fontsize',18,'Interpreter','Latex')
% grid on
% str1="/figures/figure_IRF_wefare_tau0.pdf";
% str_save1=append(root,str1);
% settt
% print(gcf,'-dpdf',str_save1)
% 
% figure(8)
% plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'--r','LineWidth',10),hold on
% plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100*log(y_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)/y_time_s_mat(end,inde_mix,1)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100*log(y_time_s_mat(1:T_irf,inde_mix2,inde_size_benefit)/y_time_s_mat(end,inde_mix,1)),'-b','LineWidth',6),hold off
% xlabel('Year','Fontsize',18,'Interpreter','Latex')
% ylabel('Output, \% deviation from status quo','Fontsize',18,'Interpreter','Latex')
% grid on
% str1="/figures/figure_IRF_output_tau0.pdf";
% str_save1=append(root,str1);
% settt
% print(gcf,'-dpdf',str_save1)
% 
% 
% figure(9)
% plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,[zeros(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100*log(i_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)/i_time_s_mat(end,inde_mix,1)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,100*log(i_time_s_mat(1:T_irf,inde_mix2,inde_size_benefit)/i_time_s_mat(end,inde_mix,1)),'-b','LineWidth',6),hold on
% xlabel('Year','Fontsize',18,'Interpreter','Latex')
% ylabel('Investment, \% deviation from status quo','Fontsize',18,'Interpreter','Latex')
% grid on
% str1="/figures/figure_IRF_investment_tau0.pdf";
% str_save1=append(root,str1);
% settt
% print(gcf,'-dpdf',str_save1)
% 
% figure(10)
% plot(-Tm1+dt:dt:0,[ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,[ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,(l_time_s_mat(1:T_irf,inde_mix,inde_size_benefit)),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf-dt,(l_time_s_mat(1:T_irf,inde_mix2,inde_size_benefit)),'-b','LineWidth',6),hold on
% xlabel('Year','Fontsize',18,'Interpreter','Latex')
% ylabel('Workers','Fontsize',18,'Interpreter','Latex')
% grid on
% str1="/figures/figure_IRF_labor_tau0.pdf";
% str_save1=append(root,str1);
% settt
% print(gcf,'-dpdf',str_save1)
% 
% 
Tm1=3;
figure(11)
plot(-Tm1+dt:dt:0,100*shift_exit_s*[exit_ss1_s*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
plot(-Tm1+dt:dt:0,100*shift_exit_n*[exit_ss1_s*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_exit_s*([exit_ss1_s+exit_0_s(inde_mix,inde_size_benefit),exit_t_s(1:T_irf/dt,inde_mix,inde_size_benefit)']/dt),'-b','LineWidth',6),hold on
plot(0:dt:T_irf,100*shift_exit_s*([exit_ss1_s+exit_0_s(inde_mix2,inde_size_benefit),exit_t_s(1:T_irf/dt,inde_mix2,inde_size_benefit)']/dt),'--r','LineWidth',6),hold on
xlabel('Year','Fontsize',18,'Interpreter','Latex')
ylabel('Exit rate, in \%','Fontsize',18,'Interpreter','Latex')
le=legend('$\vartheta+\tau=0.35$, $\tilde\tau=0.00$','$\vartheta+\tau=0.35$, $\tilde\tau=0.75$');
set(le,'Interpreter','Latex','Fontsize',16,'Location','NorthEast');
str1="/figures/figure_IRF_exit_tau0.pdf";
str_save1=append(root,str1);
grid on
settt
print(gcf,'-dpdf',str_save1)

% 
% figure(12)
% plot(-Tm1+dt:dt:0,100*shift_entry_s*[entry_ss1_s*ones(1,Tm1/dt)]/dt,'--r','LineWidth',6),hold on
% plot(-Tm1+dt:dt:0,100*shift_entry_n*[entry_ss1_s*ones(1,Tm1/dt)]/dt,'-b','LineWidth',6),hold on
% plot(0:dt:T_irf,100*shift_entry_s*([entry_t_s(1:(T_irf+dt)/dt,inde_mix,inde_size_benefit)']/dt),'--r','LineWidth',6),hold on
% plot(0:dt:T_irf,100*shift_entry_n*[entry_t_s(1:(T_irf+dt)/dt,inde_mix2,inde_size_benefit)']/dt,'-b','LineWidth',6),hold off
% xlabel('Year','Fontsize',18,'Interpreter','Latex')
% ylabel('Entry rate, in \%','Fontsize',18,'Interpreter','Latex')
% str1="/figures/figure_IRF_entry_tau0.pdf";
% str_save1=append(root,str1);
% grid on
% settt
% print(gcf,'-dpdf',str_save1)
% 
% 
% 
